{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Markov matrices\n",
    "\n",
    "A matrix $A$ is a **Markov matrix** if\n",
    "\n",
    "* Its entries are all $\\ge 0$\n",
    "* Each **column**'s entries **sum to 1**\n",
    "\n",
    "Typicaly, a Markov matrix's entries represent **transition probabilities** from one state to another.\n",
    "\n",
    "For example, consider the $2 \\times 2$ Markov matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " 0.9  0.2\n",
       " 0.1  0.8"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [0.9 0.2\n",
    "     0.1 0.8]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us suppose that this represents the fraction of people switching majors each year between math and English literature.\n",
    "\n",
    "Let\n",
    "$$\n",
    "x = \\begin{pmatrix} m \\\\ e \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "represent the number of math majors $m$ and English majors $e$.  Suppose that each year, 10% of math majors and 20% of English majors switch majors.  After one year, the new number of math and English majors is:\n",
    "\n",
    "$$\n",
    "m' = 0.9 m + 0.2 e \\\\\n",
    "e' = 0.1 m + 0.8 e\n",
    "$$\n",
    "\n",
    "But this is equivalent to a matrix multiplication!  i.e. the numbers $x'$ of majors after one year is\n",
    "\n",
    "$$\n",
    "x' = A x \\,\n",
    "$$\n",
    "\n",
    "Note that the two Markov properties are critical: we never have negative numbers of majors (or negative probabilities), and the probabilities must sum to 1 (the net number of majors is not changing: we're not including new students or people that graduate in this silly model)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eigenvalues of Markov matrices\n",
    "\n",
    "There are two key questions about Markov matrices that can be answered by analysis of their eigenvalues:\n",
    "\n",
    "* Is there a **steady state**?\n",
    "  - i.e. is there an $x_0 \\ne 0$ such that $A x_0 = x_0$?\n",
    "  - i.e. is there $\\lambda_0 = 1$ eigenvector $x_0$?\n",
    "\n",
    "* Does the system **tend toward a steady state?**\n",
    "  - i.e. does $A^n x \\to \\mbox{multiple of } x_0$ as $n \\to \\infty$?\n",
    "  - i.e. is $\\lambda = 1$ the **largest** $|\\lambda|$?\n",
    "  \n",
    "The answers are **YES** for **any Markov** matrix $A$, and **YES** for any *positive* Markov matrix (Markov matrices with entries $> 0$, not just $\\ge 0$).  For *any* Markov matrix, all of the λ satisfy $|\\lambda| \\le 1$, but if there are zero entries in the matrix we *may* have multiple $|\\lambda|=1$ eigenvalues (though this doesn't happen often in practical Markov problems)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 0.7\n",
       " 1.0"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigvals(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try just multipling it many times by a \"random\" vector and see whether it is converging to a steady state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 14.000000000000089\n",
       "  7.000000000000044"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A^100 * [17, 4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yes, it seems to be giving a vector that is not changing, which shoud be a multiple $c x_0$ of a steady-state eigenvector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 14.000000000000874\n",
       "  7.000000000000437"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cx₀ = A^1000 * [17, 4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check that this is an eigenvector of $A$ with eigenvalue $\\lambda=1$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 14.000000000000874\n",
       "  7.000000000000437"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A * cx₀"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see why, the key idea is to write the columns-sum-to-one property of Markov matrices in linear-algebra terms.  It is equivalent to the statement:\n",
    "\n",
    "$$\n",
    "\\underbrace{\\begin{pmatrix} 1 & 1 & \\cdots & 1 & 1 \\end{pmatrix}}_{o^T} A = o^T\n",
    "$$\n",
    "\n",
    "since this is just the operation that sums all of the rows of $A$.  Equivalently, if we transpose both sides:\n",
    "\n",
    "$$\n",
    "A^T o = o\n",
    "$$\n",
    "\n",
    "i.e. $o$ is an eigenvector of $A^T$ (called a **left eigenvector of A**) with eigenvalue $\\lambda = 1$.\n",
    "\n",
    "But since $A$ and $A^T$ have the **same eigenvalues** (they have the same characteristic polynomial $\\det (A - \\lambda I) = \\det (A^T - \\lambda I)$ because transposed don't change determinants), this means that $A$ **also has an eigenvalue 1** but with a **different eigenvector**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×2 adjoint(::Vector{Float64}) with eltype Float64:\n",
       " 1.0  1.0"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "o = [1,1]\n",
    "o'A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 1.0\n",
       " 1.0"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A' * o"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An eigenvector of $A$ with eigenvalue $1$ must be a basis for $N(A - I)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " -0.1   0.2\n",
       "  0.1  -0.2"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A - 1*I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By inspection, $A - I$ is singular here: the second column is -2 times the first.  So, $x_0 = (2,1)$ is a basis for its nullspace, and is a steady state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 5.551115123125783e-17\n",
       " 5.551115123125783e-17"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A - I) * [2,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check if some arbitrary starting vector $(3,0)$ tends towards this steady state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\n",
       "\\begin{pmatrix} 3.0\\\\0.0\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.7\\\\0.3\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.49\\\\0.51\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.343\\\\0.657\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.24\\\\0.76\\end{pmatrix} \\longrightarrow \\\\\n",
       "\\begin{pmatrix} 2.24\\\\0.76\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.168\\\\0.832\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.118\\\\0.882\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.082\\\\0.918\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.058\\\\0.942\\end{pmatrix} \\longrightarrow \\\\\n",
       "\\begin{pmatrix} 2.058\\\\0.942\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.04\\\\0.96\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.028\\\\0.972\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.02\\\\0.98\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.014\\\\0.986\\end{pmatrix} \\longrightarrow \\\\\n",
       "\\begin{pmatrix} 2.014\\\\0.986\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.01\\\\0.99\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.007\\\\0.993\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.005\\\\0.995\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.003\\\\0.997\\end{pmatrix} \\longrightarrow \\\\\n",
       "\\begin{pmatrix} 2.003\\\\0.997\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.002\\\\0.998\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.002\\\\0.998\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.001\\\\0.999\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.001\\\\0.999\\end{pmatrix} \\longrightarrow \\\\\n",
       "\\begin{pmatrix} 2.001\\\\0.999\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.001\\\\0.999\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.0\\\\1.0\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.0\\\\1.0\\end{pmatrix} \\longrightarrow \\begin{pmatrix} 2.0\\\\1.0\\end{pmatrix} \\longrightarrow \\\\\n",
       "$\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The following code prints a sequence of Aⁿx values\n",
    "# for n=0,1,2,… nicely formatted with LaTeX.\n",
    "\n",
    "x = [3, 0]\n",
    "pmatrix(x) = string(\"\\\\begin{pmatrix} \", round(x[1],digits=3), \"\\\\\\\\\", round(x[2],digits=3), \"\\\\end{pmatrix}\")\n",
    "buf = IOBuffer()\n",
    "println(buf, \"\\$\")\n",
    "for k = 1:6\n",
    "    print(buf, pmatrix(x), \" \\\\longrightarrow \")\n",
    "    for i = 1:4\n",
    "        x = A*x\n",
    "        print(buf, pmatrix(x), \" \\\\longrightarrow \")\n",
    "    end\n",
    "    println(buf, \"\\\\\\\\\")\n",
    "end\n",
    "println(buf, \"\\$\")\n",
    "display(\"text/latex\", String(take!(buf)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yes!  In fact, it tends to exactly $(2,1)$, because the other eigenvalue is $< 1$ (and hence that eigenvector component decays exponentially fast).\n",
    "\n",
    "An interesting property is that the **sum of the vector components is conserved** when we multiply by a Markov matrix.  Given a vector $x$, $o^T x$ is the sum of its components.  But $o^T A = o^T$, so:\n",
    "\n",
    "$$\n",
    "o^T A x = o^T x = o^T A^n x\n",
    "$$\n",
    "\n",
    "for any $n$!  This is why $(3,0)$ must tend to $(2,1)$, and not to any other multiple of $(2,1)$, because both of them sum to 3.  (The \"number of majors\" is conserved in this problem.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why no eigenvalues > 1?\n",
    "\n",
    "Why are all $|\\lambda| \\le 1$ for a Markov matrix?\n",
    "\n",
    "The key fact is that the **product AB of two Markov matrices A and B is also Markov**.  Reasons:\n",
    "\n",
    "* If $A$ and $B$ have nonnegative entries, $AB$ does as well: matrix multiplication uses only $\\times$ and $+$, and can't introduce a minus sign.\n",
    "\n",
    "* If $o^T A = o^T$ and $o^T B = o^T$ (both have columns summing to 1), then $o^T AB = o^T B = o^T$: the columns of $AB$ sum to 1.\n",
    "\n",
    "For example, $A^n$ is a Markov matrix for any $n$ if $A$ is Markov.\n",
    "\n",
    "Now, if there were an eigenvalue $|\\lambda| > 1$, the matrix $A^n$ would have to *blow up exponentially* as $n\\to \\infty$ (since the matrix times that eigenvector, or any vector with a nonzero component of that eigenvector, would blow up).  But since $A^n$ is Markov, all of its entries must be between 0 and 1.  It can't blow up!  So we must have all $|\\lambda| \\le 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " 0.666667  0.666667\n",
       " 0.333333  0.333333"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A^100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(In fact, $A^n$ is pretty boring for large $n$: it just takes in any vector and redistributes it to the steady state.)\n",
    "\n",
    "Another way of thinking about $A^{100}$ is\n",
    "$$\n",
    "A^{100} = A^{100} \\begin{pmatrix} 1 & 0 \\\\ 0 & 1 \\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    "   A^{100} \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} &\n",
    "   A^{100} \\begin{pmatrix} 0 \\\\ 1 \\end{pmatrix}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "i.e. it multiplies $A^{100}$ by each column of the identity matrix (= different possible \"starting populations\").   Because of this, each column of $A^{100}$ tends towards an eigenvector with the biggest $|\\lambda|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Can there be more than one steady state?\n",
    "\n",
    "We have just showed that we have *at least one* eigenvalue $\\lambda = 1$, and that *all* eigenvalues satisfy $|\\lambda| \\le 1$.  But can there be *more than one* independent eigenvector with $\\lambda = 1$?\n",
    "\n",
    "**Yes!** For example, the **identity matrix** $I$ is a Markov matrix, and *all* of its eigenvectors have eigenvalue $1$.  Since $Ix = x$ for *any* $x$, *every vector is a steady state* for $I$!\n",
    "\n",
    "But this does not usually happen for *interesting* Markov matrices coming from real problems.  In fact, there is a theorem:\n",
    "\n",
    "* If all the entries of a Markov matrix are $> 0$ (not just $\\ge 0$), then *exactly one* of its eigenvalues $\\lambda = 1$ (that eigenvalue has \"multiplicity 1\": $N(A-I)$ is one-dimensional), and **all other eigenvalues have** $|\\lambda| < 1$.  There is a *unique steady state* (up to an overall scale factor).\n",
    "\n",
    "I'm not going to prove this in 18.06, however."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Can the solutions oscillate?\n",
    "\n",
    "If you have a Markov matrix with zero entries, then there might be more than one eigenvalue with $|\\lambda| = 1$, but these additional solutions might be *oscillating* solutions rather than steady states.\n",
    "\n",
    "For example, consider the permutation matrix\n",
    "$$\n",
    "P = \\begin{pmatrix} 0 & 1 \\\\ 1 & 0 \\end{pmatrix}\n",
    "$$\n",
    "that simply swaps the first and second entries of any 2-component vector.\n",
    "\n",
    "If $x = (1,0)$, then $P^n x$ will oscillate forever, never reaching a steady state!  It simply oscillates between $(1,0)$ (for even $n$) and $(0,1)$ (for odd $n$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Int64}:\n",
       " 0  1\n",
       " 1  0"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = [0 1\n",
    "     1 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Vector{Vector{Int64}}:\n",
       " [1, 0]\n",
       " [0, 1]\n",
       " [1, 0]\n",
       " [0, 1]\n",
       " [1, 0]\n",
       " [0, 1]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[P^n * [1,0] for n = 0:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But this is a Markov matrix, so all $|\\lambda|$ are $\\le 1$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " -1.0\n",
       "  1.0"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigvals(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem is that the $\\lambda = -1$ eigenvalue corresponds to an oscillating solution:\n",
    "\n",
    "$$\n",
    "P^n \\begin{pmatrix} 1 \\\\ -1 \\end{pmatrix} = (-1)^n \\begin{pmatrix} 1 \\\\ -1 \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "for the eigenvector $(1,-1)$.\n",
    "\n",
    "The steady state still exists, corresponding to the eigenvector $(1,1)$:\n",
    "\n",
    "$$\n",
    "P^n \\begin{pmatrix} 1 \\\\ 1 \\end{pmatrix} =  \\begin{pmatrix} 1 \\\\ 1 \\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " -0.707107  0.707107\n",
       "  0.707107  0.707107"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = eigvecs(P) # the eigenvectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       "  1.0  1.0\n",
       " -1.0  1.0"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X ./ X[1,:]' # normalize the first row to be 1, to resemble our hand solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since $(1,0) = [(1,1) + (1,-1)]/2$, we have:\n",
    "\n",
    "$$\n",
    "P^n \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} = \\frac{1}{2} \\left[ \\begin{pmatrix} 1 \\\\ 1 \\end{pmatrix} + \n",
    "(-1)^n \\begin{pmatrix} 1 \\\\ -1 \\end{pmatrix} \\right]\n",
    "$$\n",
    "\n",
    "which alternates between $(1,0)$ and $(0,1)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another example\n",
    "\n",
    "Let's generate a random 5x5 Markov matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Matrix{Float64}:\n",
       " 0.410618  0.306837  0.410031  0.707623   0.290909\n",
       " 0.307687  0.22414   0.676996  0.0455438  0.904309\n",
       " 0.999213  0.714056  0.357485  0.913338   0.715352\n",
       " 0.647026  0.995701  0.789245  0.577309   0.391341\n",
       " 0.73899   0.967951  0.914835  0.565266   0.447786"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = rand(5,5) # random entries in [0,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×5 Matrix{Float64}:\n",
       " 3.10353  3.20869  3.14859  2.80908  2.7497"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(M,dims=1) # not Markov yet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Matrix{Float64}:\n",
       " 0.132307   0.0956271  0.130227  0.251906   0.105797\n",
       " 0.0991408  0.0698541  0.215016  0.0162131  0.328876\n",
       " 0.32196    0.222539   0.113538  0.325138   0.260157\n",
       " 0.20848    0.310314   0.250666  0.205515   0.142322\n",
       " 0.238113   0.301666   0.290554  0.201228   0.162849"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = M ./ sum(M,dims=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×5 Matrix{Float64}:\n",
       " 1.0  1.0  1.0  1.0  1.0"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(M,dims=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{ComplexF64}:\n",
       "   -0.17776058595953462 + 0.0im\n",
       "    -0.1352931760939033 + 0.0im\n",
       " -0.0014416561523480091 - 0.07745841517651891im\n",
       " -0.0014416561523480091 + 0.07745841517651891im\n",
       "     1.0000000000000004 + 0.0im"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigvals(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{Float64}:\n",
       " 0.17776058595953462\n",
       " 0.1352931760939033\n",
       " 0.07747183006822272\n",
       " 0.07747183006822272\n",
       " 1.0000000000000004"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "abs.(eigvals(M))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{Float64}:\n",
       " 0.05662283043686728\n",
       " 0.13759834468209436\n",
       " 0.3955727782747076\n",
       " 0.09136601599437516\n",
       " 0.3188400306119557"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = rand(5)\n",
    "x = x / sum(x) # normalize x to have sum = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{Float64}:\n",
       " 0.14590618751218656\n",
       " 0.15842031877820567\n",
       " 0.24194675375641556\n",
       " 0.2186167846876927\n",
       " 0.23510995526549555"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M^100 * x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.0, 0.9999999999999962)"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(x), sum(M^100 * x) # still = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{ComplexF64}:\n",
       "  0.1459061875121874 + 0.0im\n",
       " 0.15842031877820623 + 0.0im\n",
       " 0.24194675375641644 + 0.0im\n",
       " 0.21861678468769355 + 0.0im\n",
       " 0.23510995526549633 + 0.0im"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "λ, X = eigen(M)\n",
    "X[:,end] / sum(X[:,end]) # eigenvector for λ=1, normalized to sum=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, $M^n x$ is approaching a steady-state ($\\lambda = 1$) eigenvector of $M$ as $n$ grows large."
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.2"
  },
  "widgets": {
   "state": {
    "e53e5f7b-c65e-4676-a564-3f8ee40c11c0": {
     "views": [
      {
       "cell_index": 13
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
